#How Do Observers Assess Resolve?
#July 7, 2018
#BJPS Replication 7.R

#This file contains the replication code needed to replicate the analysis from Supplementary Appendix §1.42 (Average marginal treatment interaction effects)

#To generate the necessary dataframes, either run BJPS Replication 1.R, or load the following saved object:
library(here) #Avoids needing to setwd()
dat3 <- get(load(file="Resolve Conjoint Cleaned No US 070718.RData"))

#Load FindIt package
library(FindIt)

### Note: lines 14-34 take a while, even on cluster computing, so we've also saved the model objects directly; the models below were generated on the Harvard-MIT Data Center RCE

#### Estimate model with up to 3-way interactions between factors
conjoint.mod <- FindIt(model.treat=Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, nway=3, data=dat, type="binary", treat.type="multiple")
print("Model successfully estimated! Rejoice!")
#Save model
save(conjoint.mod, file="Resolve Conjoint SVM Mod.RData")

#Now, generate full factorial design matrix
full.FindIt <- function(object){
    ## Make the full factorial design matrix as the target population.
    comb.list <- apply(object$treat.orig,2,unique)
    full <- as.data.frame(expand.grid(comb.list))
    invisible(full)
}

full <- full.FindIt(conjoint.mod) #Produces full factorial design matrix
pred.dat <- predict(conjoint.mod, newdata=full, sort=FALSE)$data #Predicted potential outcomes
print("Predicted potential outcomes complete! Hurray!")
#Save model
save(pred.dat, file="Resolve Conjoint SVM Pred Dat.RData")

## Load data in case you don't want to run lines 14-34
conjoint.mod <- get(load("Resolve Conjoint SVM Mod.RData"))
pred.dat <- get(load("Resolve Conjoint SVM Pred Dat.RData"))

##The INT function. This is taken from the FindIt package version 0.6, by Naoki Egami, Marc Ratkovic and Kosuke Imai. Newer versions of the package use a slightly different estimation methods.
INT <- function (object, target.data, column, dist = "target", base, 
    sort = TRUE, compare = FALSE, order = 2) 
{
    if (all(class(object) != "FindIt")) {
        warning("the class of object needs to be FindIt.")
    }
    Restriction <- NULL
    measure <- "TIE"
    if (compare == TRUE) {
        if (missing(order)) {
            warning("Need to specify the order")
        }
    }
    coefs <- object$coefs.orig
    names(coefs) <- c("Intercept", object$name.t)
    outcome.type <- object$type
    if (dist == "unique") {
        data <- cbind(object$y.orig, object$treat.orig)
        colnames(data)[-1] <- colnames(object$treat.orig)
        data <- unique(data, MARGIN = 1)
    }
    if (dist == "sample") {
        data <- cbind(object$y.orig, object$treat.orig)
        colnames(data)[-1] <- colnames(object$treat.orig)
    }
    if (dist == "target") {
        data <- target.data
        print("Using Data representing the target population.")
    }
    else {
        print("Using the Data used to fit the model.")
    }
    if (!missing(column) & !missing(base)) {
        for (i in 1:length(column)) {
            ind.column <- which(colnames(data) == column[i])
            data[, ind.column] <- relevel(data[, ind.column], 
                base[i])
        }
    }
    if (!missing(column)) {
        ind <- c()
        for (i in 1:length(column)) {
            ind[i] <- which(colnames(data) == column[i])
        }
        column <- column[order(ind, decreasing = FALSE)]
    }
    base <- 1
    Marginal <- function(data, base, Restriction) {
        Effect.list <- list()
        Name.list <- list()
        range.mar <- c()
        for (i in 1:(ncol(data) - 1)) {
            columnM <- colnames(data)[i + 1]
            A <- tapply(data[, 1], data[, columnM], mean, simplify = FALSE)
            Effect1 <- unlist(A)
            if (base == "min") {
                base <- which(Effect1 == min(Effect1))
            }
            TEffect <- Effect1 - Effect1[base]
            range.mar[i] <- max(TEffect) - min(TEffect)
            Effect.list[[i]] <- TEffect
            Name.list[[i]] <- paste(colnames(data)[(i + 1)], 
                names(TEffect), sep = "_")
        }
        names(range.mar) <- colnames(data)[2:ncol(data)]
        name <- unlist(Name.list)
        Treatment.Effect1 <- unlist(Effect.list)
        Treatment.Effect1 <- as.data.frame(Treatment.Effect1)
        colnames(Treatment.Effect1) <- "AMTE"
        rownames(Treatment.Effect1) <- name
        output <- list(Treatment.Effect1 = Treatment.Effect1, 
            range = range.mar)
        return(output)
    }
    if (compare == FALSE & missing(column)) {
        Treatment.Effect1 <- Marginal(data = data, base = base, 
            Restriction = Restriction)$Treatment.Effect1
        Treatment.Effect <- Treatment.Effect1
    }
    if (compare == TRUE & order == 1) {
        Marginal.range <- Marginal(data = data, base = base, 
            Restriction = Restriction)$range
        print("Range of Marginal Effects")
        return(Marginal.range)
    }
    comb.prem <- function(data, column, Restriction, order) {
        A <- tapply(data[, 1], data[, column], mean, simplify = FALSE)
        A2 <- tapply(data[, 1], data[, column], mean, simplify = TRUE)
        A[is.na(A2)] <- NA
        Effect1 <- unlist(A)
        if (base == "min") {
            base <- which(Effect1 == min(Effect1))
        }
        Comb.Effect <- Effect1 - Effect1[base]
        Treatment.Effect <- cbind(Comb.Effect, expand.grid(dimnames(A)))
        Treatment.Effect <- na.omit(Treatment.Effect)
        Comb.Effect <- Treatment.Effect[, 1]
        Combination.name <- Treatment.Effect[, -1]
        Treatment.Effect.print <- Treatment.Effect
        for (j in 2:ncol(Treatment.Effect)) {
            Treatment.Effect[, j] <- paste(colnames(Treatment.Effect)[j], 
                Treatment.Effect[, j], sep = "_")
        }
        Treatment.Effect1 <- Marginal(data = data, base = base, 
            Restriction = Restriction)$Treatment.Effect1
        Main.comb <- c()
        for (i in 1:nrow(Treatment.Effect)) {
            main <- c()
            for (j in 1:(ncol(Treatment.Effect) - 1)) {
                main[j] <- Treatment.Effect1[rownames(Treatment.Effect1) == 
                  Treatment.Effect[i, (j + 1)], 1]
            }
            Main.comb[i] <- sum(main)
        }
        Sum.Mar.Effect <- Main.comb - Main.comb[base]
        if (order == 2) {
            TIE <- round(Comb.Effect - Sum.Mar.Effect, digits = 8)
        }
        if (order == 3) {
            data.column.three <- data[, column]
            Two.way.comb.list <- list()
            for (j in 1:3) {
                column.com.three <- c()
                Two.way.comb <- list()
                column.com.three <- colnames(data.column.three)[c(combn(length(data.column.three), 
                  2)[, j])]
                A.Three <- tapply(data[, 1], data[, column.com.three], 
                  mean, simplify = FALSE)
                A2.Three <- tapply(data[, 1], data[, column.com.three], 
                  mean, simplify = TRUE)
                A.Three[is.na(A2.Three)] <- NA
                EffectTwo <- unlist(A.Three)
                Two.way.comb <- EffectTwo - EffectTwo[base]
                Two.way.comb2 <- cbind(Two.way.comb, expand.grid(dimnames(A.Three)))
                for (k in 2:3) {
                  Two.way.comb2[, k] <- paste(colnames(Two.way.comb2)[k], 
                    Two.way.comb2[, k], sep = "_")
                }
                colnames(Two.way.comb2) <- c("Two.way.comb", 
                  "first", "second")
                Two.way.comb.list[[j]] <- na.omit(Two.way.comb2)
            }
            Two.way.comb.Final <- do.call(rbind, Two.way.comb.list)
            Two.Way.sum <- c()
            for (i in 1:nrow(Treatment.Effect)) {
                two.sum <- c()
                for (j in 1:nrow(Two.way.comb.Final)) {
                  ind.two <- sum(is.element(Two.way.comb.Final[j, 
                    2:3], Treatment.Effect[i, 2:4]))
                  if (ind.two == 2) {
                    two.sum[j] <- Two.way.comb.Final[j, 1]
                  }
                  else {
                    two.sum[j] <- 0
                  }
                }
                Two.Way.sum[i] <- sum(two.sum)
            }
            Sum.Two.Effect <- Two.Way.sum - Two.Way.sum[base]
            TIE <- round(Comb.Effect - Sum.Two.Effect + Sum.Mar.Effect, 
                digits = 8)
        }
        range.change <- max(TIE) - min(TIE)
        Order.data <- cbind(Sum.Mar.Effect, Comb.Effect)
        Order.data <- Order.data[order(Order.data[, 1], decreasing = TRUE), 
            ]
        Order.data <- as.data.frame(Order.data)
        Order.data$order.main <- seq(1:nrow(Order.data))
        Order.comb <- Order.data[order(Order.data[, 2], decreasing = TRUE), 
            "order.main"]
        Order.diff <- Order.comb - seq(1:nrow(Order.data))
        max.order.diff <- max(abs(Order.diff))
        return(list(A = A, Treatment.Effect = Treatment.Effect, 
            Comb.Effect = Comb.Effect, Combination.name = Combination.name, 
            Sum.Mar.Effect = Sum.Mar.Effect, TIE = TIE, range.change = range.change, 
            max.order.diff = max.order.diff))
    }
    if (compare == FALSE) {
        if (!missing(column)) {
            X <- comb.prem(data = data, column = column, Restriction = Restriction, 
                order = order)
            A <- X$A
            Treatment.Effect <- X$Treatment.Effect
            Comb.Effect <- X$Comb.Effect
            Combination.name <- X$Combination.name
            Sum.Mar.Effect <- X$Sum.Mar.Effect
            TIE <- X$TIE
            range.change <- X$range.change
        }
    }
    else {
        data.column <- data[, -1]
        column.com <- list()
        Range.com <- list()
        for (j in 1:choose(length(data.column), order)) {
            column.com[[j]] <- colnames(data.column)[c(combn(length(data.column), 
                order)[, j])]
            if (measure == "TIE") {
                Range.com[[j]] <- comb.prem(data = data, column = column.com[[j]], 
                  Restriction = Restriction, order = order)$range.change
            }
            if (measure == "Order.Diff") {
                Range.com[[j]] <- comb.prem(data = data, column = column.com[[j]], 
                  Restriction = Restriction, order = order)$max.order.diff
            }
        }
        column.com <- as.data.frame(matrix(unlist(column.com), 
            ncol = order, byrow = TRUE))
        Compare <- as.data.frame(cbind(column.com, unlist(Range.com)))
        Compare <- Compare[order(Compare[, (order + 1)], decreasing = TRUE), 
            ]
        if (measure == "TIE") {
            colnames(Compare) <- c(colnames(Compare)[1:order], 
                "range.TIE")
        }
        if (measure == "Order.Diff") {
            colnames(Compare) <- c(colnames(Compare)[1:order], 
                "max.Order.diff")
        }
        return(Compare)
    }
    if (!missing(column)) {
        if (length(column) == 2) {
            Treatment.coefs.name <- Treatment.Effect[, 2]
            if (ncol(Treatment.Effect) >= 3) {
                for (i in 1:nrow(Treatment.Effect)) {
                  for (j in 3:ncol(Treatment.Effect)) {
                    Treatment.coefs.name[i] <- paste(Treatment.coefs.name[i], 
                      Treatment.Effect[i, j], sep = ".")
                  }
                }
            }
            Treatment.coefs.name <- gsub(" ", ".", Treatment.coefs.name)
        }
    }
    if (missing(column)) {
        Treatment.coefs.name <- rownames(Treatment.Effect1)
        Treatment.coefs.name <- gsub(" ", ".", Treatment.coefs.name)
    }
    if (missing(column)) {
        Coefficients <- c()
        for (i in 1:length(Treatment.coefs.name)) {
            if (Treatment.coefs.name[i] %in% names(coefs)) {
                Coefficients[i] <- coefs[names(coefs) == Treatment.coefs.name[i]]
            }
            else {
                Coefficients[i] <- NA
            }
        }
        if (outcome.type == "binary") {
            Coefficients <- Coefficients/2
        }
    }
    if (!missing(column)) {
        if (length(column) == 2) {
            Coefficients <- c()
            for (i in 1:length(Treatment.coefs.name)) {
                if (Treatment.coefs.name[i] %in% names(coefs)) {
                  Coefficients[i] <- coefs[names(coefs) == Treatment.coefs.name[i]]
                }
                else {
                  Coefficients[i] <- NA
                }
            }
            if (outcome.type == "binary") {
                Coefficients <- Coefficients/2
            }
        }
    }
    if (!missing(column)) {
        if (length(column) == 2) {
            Treatment.Effect.print <- cbind(Comb.Effect, TIE, 
                Combination.name)
            colnames(Treatment.Effect.print) <- c("ATCE", "AMTIE", 
                colnames(Combination.name))
            Main.matrix <- cbind(Sum.Mar.Effect, TIE, Combination.name)
            colnames(Main.matrix) <- c("Sum of AMTEs", "AMTIE", 
                colnames(Combination.name))
            Inequality.Cor <- cor(Sum.Mar.Effect, TIE)
            Relative.change <- cbind(TIE, Combination.name)
            colnames(Relative.change) <- c("AMTIE", colnames(Combination.name))
            Relative.change <- Relative.change[order(Relative.change[, 
                1], decreasing = TRUE), ]
        }
        if (length(column) >= 3) {
            Treatment.Effect.print <- cbind(Comb.Effect, TIE, 
                Combination.name)
            colnames(Treatment.Effect.print) <- c("ATCE", "AMTIE", 
                colnames(Combination.name))
            Main.matrix <- cbind(Sum.Mar.Effect, TIE, Combination.name)
            colnames(Main.matrix) <- c("Sum of AMTEs", "AMTIE", 
                colnames(Combination.name))
            Relative.change <- cbind(TIE, Combination.name)
            colnames(Relative.change) <- c("AMTIE", colnames(Combination.name))
            Relative.change <- Relative.change[order(Relative.change[, 
                1], decreasing = TRUE), ]
        }
    }
    else {
        Treatment.Effect.print <- Treatment.Effect1
        Man.matrix <- NULL
    }
    Treatment.Effect.print <- as.data.frame(Treatment.Effect.print)
    if (!missing(column)) {
        Main.matrix <- as.data.frame(Main.matrix)
    }
    if (sort == TRUE) {
        Treatment.Effect.print1 <- Treatment.Effect.print
        Treatment.Effect.print <- Treatment.Effect.print[order(Treatment.Effect.print[, 
            1], decreasing = TRUE), ]
        if (!missing(column)) {
            Main.matrix <- Main.matrix[order(Main.matrix$Sum, 
                decreasing = TRUE), ]
            Treatment.Effect.print <- cbind(Treatment.Effect.print[, 
                1:2], Treatment.Effect.print[, 3:ncol(Treatment.Effect.print)])
        }
    }
    if (!missing(column)) {
        if (sort) {
            return(list(`Range of AMTIE` = range.change, AMTIE = Relative.change, 
                ATCE = Treatment.Effect.print, `Sum of AMTEs` = Main.matrix))
        }
        else {
            return(list(`Range of AMTIE` = range.change, AMTIE = Relative.change, 
                ATCE = Treatment.Effect.print, `Sum of AMTEs` = Main.matrix))
        }
    }
    else {
        return(list(Treatment.Effect = Treatment.Effect.print))
    }
}

#Now, calculate AMTIE ranges:

compare1 <- INT(conjoint.mod, target.data=pred.dat, compare=TRUE, order=1) #Main effects
compare2 <- INT(conjoint.mod, target.data=pred.dat, compare=TRUE, order=2) #Two-way interactions
compare3 <- INT(conjoint.mod, target.data=pred.dat, compare=TRUE, order=3) #Three-way interactions

##Calculate barplot data
sorted.compare1 <- sort(compare1, decreasing=TRUE)
compare.names1 <- names(sorted.compare1)

sorted.compare2 <- compare2[order(-compare2[,3]),] 
compare.names2 <- paste(eval(sorted.compare2[,1]), eval(sorted.compare2[,2]), sep = ":")

sorted.compare3 <- compare3[order(-compare3[,4]),] 
compare.names3 <- paste(eval(sorted.compare3[,1]), eval(sorted.compare3[,2]), eval(sorted.compare3[,3]), sep = ":")

#Adjust names
compare.names1r <- compare.names1
compare.names1r <- sub("prevActLdr", "Past Actions",compare.names1r)
compare.names1r <- sub("curBehavior", "Current Behavior",compare.names1r)
compare.names1r <- sub("milService", "Military Service",compare.names1r)
compare.names1r <- sub("RegimeType", "Regime Type",compare.names1r)
compare.names1r <- sub("actor", "Ally",compare.names1r)
compare.names1r <- sub("NewLeader", "New Leader",compare.names1r)
compare.names1r <- sub("prevOpponent", "Prev Opponent",compare.names1r)
compare.names1r <- sub("maleLeader", "Male Leader",compare.names1r)
compare.names1r <- sub("prevRole", "Prev Role",compare.names1r)

compare.names2r <- compare.names2
compare.names2r <- sub("prevActLdr", "Past Actions",compare.names2r)
compare.names2r <- sub("curBehavior", "Current Behavior",compare.names2r)
compare.names2r <- sub("milService", "Military Service",compare.names2r)
compare.names2r <- sub("RegimeType", "Regime Type",compare.names2r)
compare.names2r <- sub("actor", "Ally",compare.names2r)
compare.names2r <- sub("NewLeader", "New Leader",compare.names2r)
compare.names2r <- sub("prevOpponent", "Prev Opponent",compare.names2r)
compare.names2r <- sub("maleLeader", "Male Leader",compare.names2r)
compare.names2r <- sub("prevRole", "Prev Role",compare.names2r)

compare.names3r <- compare.names3
compare.names3r <- sub("prevActLdr", "Past Actions",compare.names3r)
compare.names3r <- sub("curBehavior", "Current Behavior",compare.names3r)
compare.names3r <- sub("milService", "Military Service",compare.names3r)
compare.names3r <- sub("RegimeType", "Regime Type",compare.names3r)
compare.names3r <- sub("actor", "Ally",compare.names3r)
compare.names3r <- sub("NewLeader", "New Leader",compare.names3r)
compare.names3r <- sub("prevOpponent", "Prev Opponent",compare.names3r)
compare.names3r <- sub("maleLeader", "Male Leader",compare.names3r)
compare.names3r <- sub("prevRole", "Prev Role",compare.names3r)

#### Figure 12a: One-way AMTIEs
compare.barplot.oneway <- c(as.data.frame(sorted.compare1)$sorted.compare1)
compare.barplot.reversed.oneway <- rev(compare.barplot.oneway)

barplot.names.oneway <- compare.names1r
barplot.names.reversed.oneway <- rev(barplot.names.oneway)
par(mar = c(4,12,2,2))
pos <- barplot(height=compare.barplot.reversed.oneway, horiz=TRUE, xlim=c(0,.2), names.arg = barplot.names.reversed.oneway, las=2, tck=0.01, col="black", xaxt="n", main= "One-way effects", xlab="Ranges of the one-Way AMTIE", family="sans")
axis(side=2, at=pos,labels=FALSE, tck=-0.02) ##adding tick marks next to the labels
axis(side=1) ##adding labels to the x axis 
box()
abline(v=seq(0,0.2,by=0.05), lty=3, col=rgb(0,0,0,0.5))

##### Figure 12b: Two-way AMTIEs
length(which(sorted.compare2[,3]>0)) #25 are nonzero
length(which(sorted.compare2[,3]==0)) #30 are zero

compare.barplot.twoway <- sorted.compare2[, 3][which(sorted.compare2[,3]>0)]
compare.barplot.reversed.twoway <- rev(compare.barplot.twoway)

barplot.names.twoway <- compare.names2r[which(sorted.compare2[,3]>0)]
barplot.names.reversed.twoway <- rev(barplot.names.twoway)

par(mar = c(4,12,2,2))
pos <- barplot(height=compare.barplot.reversed.twoway, horiz=TRUE, xlim=c(0,.2), names.arg = barplot.names.reversed.twoway, las=2, tck=0.01, col="black", xaxt="n", main= "Two-way effects", xlab="Ranges of the two-way AMTIE", family="sans")
axis(side=2, at = pos, labels=FALSE, tck=-0.02) ##adding tick marks next to the labels
axis(side=1) ##adding labels to the x axis
box()
abline(v=seq(0,0.2,by=0.05), lty=3, col=rgb(0,0,0,0.5))

###### Figure 13c: Three-way AMTIEs

length(which(sorted.compare3[,4]>0)) #5 are nonzero
length(which(sorted.compare3[,4]==0)) #160 are zero

compare.barplot.threeway <- sorted.compare3[, 4][which(sorted.compare3[,4]>0)]
compare.barplot.reversed.threeway <- rev(compare.barplot.threeway)

barplot.names.threeway <- compare.names3r[which(sorted.compare3[,4]>0)]
barplot.names.reversed.threeway <- rev(barplot.names.threeway)

par(mar = c(4,12,2,2))
pos <- barplot(height=compare.barplot.reversed.threeway, horiz=TRUE, xlim=c(0,.2), names.arg = barplot.names.reversed.threeway, las=2, tck=0.01, col="black", xaxt="n", main= "Three-way effects", xlab="Ranges of the three-way AMTIE", family="sans")
axis(side=2, at = pos, labels=FALSE, tck=-0.02) ##adding tick marks next to the labels
axis(side=1) ##adding labels to the x axis 
box()
abline(v=seq(0,0.2,by=0.05), lty=3, col=rgb(0,0,0,0.5))

####### Figure 13d: Largest AMTIES

compare.barplot <- c(0, as.data.frame(sorted.compare1)$sorted.compare1[1:5], 0, 0, sorted.compare2[1:5, 3], 0, 0, sorted.compare3[1:5, 4])
compare.barplot.reversed <- rev(compare.barplot)

barplot.names <- c("One-way effects", compare.names1r[1:5], "", "Two-way effects", compare.names2r[1:5], "",
                   "Three-way effects", compare.names3r[1:5])

barplot.names.reversed <- rev(barplot.names)

par(mar = c(4,12,2,2))
pos <- barplot(height=compare.barplot.reversed, horiz=TRUE, xlim=c(0,.2), 
               names.arg = barplot.names.reversed, las=2, tck=0.01, col="black", xaxt="n", 
               main= "Largest AMTIEs", xlab="Ranges of the k-way AMTIE", family="sans")
axis(side=2, at = pos, labels=FALSE, tck=-0.02) ##adding tick marks next to the labels
axis(side=1) ##adding labels to the x axis (need to do this here rather in the barplot call due to the label rotation occuring in "las")
box()
abline(v=seq(0,0.2,by=0.05), lty=3, col=rgb(0,0,0,0.5))
